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ABSTRACT 

Context. The physical conditions in a collapsing cloud can be traced by observations of molecular lines. To correctly interpret these 
observations the abundance distributions of the observed species need to be derived. The chemistry in a collapsing molecular cloud is 
not in a steady state as the density and temperature evolve. We therefore need to follow chemical reactions, both in the gas phase and 
on dust grains, as well as gas-grain interactions, to predict the abundance distributions. 

Aims. Our aim is to model the abundances of molecules, in the gas phase and on grain mantles in the form of ice, from prestellar 
core collapse to disk formation. We want to investigate the need for grain surface reactions and compare our results with observed 
abundances, column densities, and ice-mantle compositions. 

Methods. We use a 2-dimensional hydrodynamical simulation as a physical model from which we take the density, temperature, and 
the flow of the gas. Trace particles, moving along with the gas, are used to follow the chemistry during prestellar core collapse and 
disk formation. We calculated abundance profiles and column densities for various species. The evolution of these abundances and 
the composition of ices on grain mantles were compared to observations and we tested the influence of grain surface reactions on the 
abundances of species. We also investigated the initial abundances to be adopted in more detailed modeling of protoplanetary disks 
by following the chemical evolution of trace particles accreting onto the disk. 

Results. Fractional abundances of HCO + , N 2 H + , H 2 CO, HC3N, and CH3OH from our model with grain surface reactions provide 
a good match to observations, while abundances of CO, CS, SO, HCN, and HNC show better agreement without grain surface 
reactions. The observed mantle composition of dust grains is best reproduced when we include surface reactions. The initial chemical 
abundances to be used for detailed modeling of a protoplanetary disk are found to be different from those in dark interstellar clouds. 
Ices with a binding energy lower than about 1200 K sublimate before accreting onto the disk, while those with a higher binding energy 
do not. 

Key words. Astrochemistry - Hydrodynamics - Stars:formation - ISM:molecules - ISM clouds 
1. Introduction and N2H + deplete at significantly higher number densities of 



Star and planet formation can be studied through emission lines | Rodge rs & Charnie^l J2003h modeled chemical abundances 
of molecules in the gas phase and from absorption lines of ices fa ^ enyel surrounding new born stars gradually heated by 



g . on dust grains. In star-forming regions, the chemistry is not the embedded protostar and£^Sl d2004) followed the chem- 
■ in steady st ate as the phys ical conditions change throughout ^ evolution duri the coll of a B onnor-Ebert sphere, 
the process <Lee et alj | 2004ft The chemical processes m star- Generall th find a hase of dep ietion followed by subli- 
forming regions are still far from understood an d give rise to mation due tQ ±Q risi te ature in the envel of the 
a^e number of different molecular species (e g.lCazauxjMjdJ tostar Different ies sublimate at different radii from 
|20Q3J). In particular, our knowledge about gas-dust interactions ^ central Qstar dug tQ differences m the bindi e ies 
and the chemical reactions on grains surfaces is still very limited. of iefj ontQ dugt in mantlef , These smdies used Qne . 
However, we need to understand these processes if we want to dimensional time -evolving models for the density and temper- 
correct y asses the physical conditions in star-forming regions ^ distribution B following fluid elements or massless trace 
and follow the chemical evolution from cold molecular cores to- ides moyi al wkh the radiall in _ falli abun . 
wards stars and planetary systems. dances couM bg deriyed ag function of distance from ^ proto _ 

Previous studies of the chemical evolution during star forma- star. No surface reactions were in cluded in both of these sim- 

tion have (mai nly) focused on pre stellar core formation and con- ulations. iGarrod & Herbstl Q2006) presented a model for hot 

traction (e.g., Aikawa et al. 2005). Freeze-out of neutral species cores which included surface reactions. However, the under- 

is an important process in these dark clouds resulting in chemi- lying physical model was relatively simple: consisting of a 

cal differentiation. This has be en observed i n numerous prestel- phase of free-fall isothermal collapse and a phase of gradual 

lar co res and dark clouds, e.g. jTafalla et alj (2002): P agani et alj heating of the gas an d dust by a central Young Stellar Object 

(120051) . These authors found that CO and sulfur-bearing species (Y SO). | Aikawa et alJ (1 2008) used the same chemical model 

(e.g., CS, SO) deplete at larger radii (lower densities) while NH3 as IGarrod & Herbstl d2006l) but adopted the results of a one- 
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dimensional radiation-hydrodynamics calculation as the phys- 
ical model of the core. They found that organic species were 
mainly formed on grain surfaces at temperatures of 20 - 40 K. 
These organic species then sublimated into the gas phase as the 
gas entered the inner hot (T > 100 K) region surrounding the 
protostar. 

Separate studies have focused on the chemi cal evolution 
in sta t ic 2-d i mensional protoplanetary disks (e.g.. Aika wa et all 
119961 Il997t lAikawa & Herbstl 1200 lb . In the cooler midplane 
of the disk various molecules freeze-out again and accumulate 
in ice mantles. Due to the strong radiation field from the pro- 
tostar ion-molecule reactio ns become important in the upper 
disk layers and above (e.g.. Ivan Dishoeck & Hogerheiide 1999; 
Ivan Dishoeck et al.ll998: Jonkheid et al. 20Ql). 



J0rgensen et al.l d2002l) derived CO abundances, using Monte 



Carlo line radiative transfer modeling, for a sample of 18 low- 
mass protostars and prestellar cores. An increase in the CO abun- 
dance from Class to Class I objects was observed (classifica- 
tion a ccording to Lada 119871: lAdams et alJ 119871: [Andre et al. 
1993). Class I objects showed almost normal cosmic fractional 
CO abundances (X co = «(CO)/n(H 2 ) ~ 2xl0~ 4 , with n(CO) and 
n(H2) the number densities of CO and H2 respectively), while for 
Class o bjects the abundances were alm ost an order of magni- 
tude lower. J0rgens e~n et al J (12004. 2005a) expanded on this work 
and constrained the abundances for a range of species in the en- 
velopes of the same 18 objects. They found various trends for the 
abundances of species as function of envelope mass, reflecting 
the fact that the chemistry is not in steady state during prestellar 
core collapse. In general, it is thought that the envelopes of Class 
I objects have a lower density and higher temperature than Class 
objects. When the temperature rises in the region surrounding 
the protostar this leads to the return of species, like CO, by sub- 
limation to the gas phase. At the same time, further out in the en- 
velope the abundances also increase as the freeze-out timescale 
gets longer due to a lower density. The picture emerging from 
these and other observations is that freeze-out and desorption are 
important processes for the evolution of chemical abundances 
from Class to Class I objects. 

Molecular lines of depleted species in Class and Class I 
objects can be account ed for with "drop"-abundance profiles, 
J0rge nsen et al.l (|2005b). In the inner and outer regions of the 
envelope there is no depletion (due to a high temperature, low 
density respectively), whereas at intermediate distances there is 
a drop in the abundance due to freeze-out. HCO + shows a clear 
relationship with CO, while N2H + abundances remain almost 
constant in the envelopes. Both of these effects can be explained 
by the formation and destruction channels of these species which 
both involve CO . CH3OH abund ances for the same sample of 
18 objects from |j0rgensen et aLl require in some cases abun- 
dance jumps in the inner regions. H2CO abundances were best 
fitted with "drop"-abundance profiles, although an inner jump 
for H2CO could not be ruled out for most sources. Abundances 
of complex species like CH3OH are not directly affected by the 
longer freeze-out timescales in the envelope because of their 
high binding energies. These species can only return (in large 
amounts) to the gas phase by grain destruction mechanism (e.g., 
in shocks) or by sublimation if the temperature is high enough 
(T > 60 K, depending on the species). 

The composition of ices on grain mantles is also of interest. 
Grai n mantles provide surfaces where species can meet and react 
(e.g. jAllen & Robinson! [19771) . Abundances of certain species, 
a typical example is CH3OH, can increase dram atically when 
comp ared to their formation in the gas phase (e.g., Garr odet alJ 
|2007|) . Also, the mantle composition has an effect on the binding 



energy of species onto grains and thus couples to the gas-grain 
interactions and gas phase chemistry. Modeled ice fractions can 
be compared to the observed ice mantle compositions in order to 
constrain which surface reactions take place. 

The evolution of chemical abundances during prestellar core 
collapse may also provide us with a so called "chemical clock". 
It is still unclear wether such a clock exists. For example the 
use of su lfur-bearing spec i es as a chemical clock was investi- 
gated by Buckle & Fuller (2003), but no definitive conclusion 
has been reached. Modeled abundances can in principle also be 
used to constrain hydrodynamical models of star formation if 
the chemistry is well enough understood. For example, high gas 
phase abundances of complex species formed on grain mantles, 
require the temperature to be at least higher than the sublimation 
tempe ratures of these species. 

In iBrinch. van W eeren. & Hogerh elidd d200 8b). from now 
on called Paper I, we followed for the first time the 2- 
dimensional time evolution of the CO abundance during prestel- 
lar core collapse and disk formation. A 2-dimensional hydro- 
dynamical simulation was used to track the temperature, den- 
sity, and follow particle motions. We also modeled line profiles 
for CO and its isotopologue C ls O using the line excitation and 
radiat ive transfer code RATRAN dHogerheiide & van der Tata 
2000). In particular, the anti-correlation of the fractional CO 
abundance with envelope mass was reproduced. However, only a 
simple model for calculating the CO abundance was used. In this 
second paper we present a method for modeling the full chem- 
ical evolution of a collapsing prestellar core in more than one 
dimension. We replace the model for CO freeze-out/sublimation 
by a gas phase chemical network which also includes gas-grain 
interactions and grain surface reactions. This is necessary in or- 
der to derive the abundances of species which are more strongly 
influenced by gas phase or grain surface reactions. We apply this 
method to a collapsing one-solar mass rotating prestellar core 
with the physical conditions given by the 2-dimensional hydro- 
dynamical simulation from Paper I. 

Apart from the reasons mentioned above there are several 
important motivations for this work. First, with the upcoming 
higher resolution (< 1") observations of ALMA it becomes im- 
portant to extend the chemical modeling towards two and three 
dimensions as one-dimensional radial abundance distributions 
will no longer provide a good description of the observations. 
The resulting abundance distributions from our model can then 
be used in line excitation and radiative transfer codes to provide 
spatially resolved line profiles over the envelope and disk sur- 
rounding protostellar objects. This is needed to correctly derive 
the physical c onditions around protostellar objects, s uch as the 
velocity field dBrinch. Hogerheiide. & Richling||2008d) . 

Second, our modeling also provides important constraints 
for the initial abundances of sp ecies to be adopted i n more 
detailed disk modeling (e.g., ISemenov et a D 12004 [2006; 
Pavl vuchenkov et aT]|2007l) as in our simulation we can track the 
chemical evolution of gas and ice accreting onto the disk (see 
Section[3]l. 

The method presented in this paper for modeling the chem- 
ical abundances is generally applicable and does not depend 
on the particular hydrodynamical simulation we have used. 
Our method can also handle one-dimensional as well as 3- 
dimensional hydrodynamical simulations as long as values for 
the density, temperature, velocity field, and the radiation field 
are provided. 

The layout of this paper is as follows. In Section[2]we shortly 
summarize the hydrodynamical simulation we use as input for 
our chemistry. We also present the chemical network to com- 
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pute the gas phase abundances and ice composition. In Section[3] 
we describe our method to couple the output from the hydrody- 
namical simulation to the chemistry. In Section[4]we present our 
results and compare them with observations. This is followed by 
a discussion and conclusions in Section[5]and[6] 



2. Model 

2.1. Physical model: a 2-dimensional hydrodynamical 
simulation 

We use the same physical model of a collapsing prestellar core 
as in Pape r I. The grid-based 2-d i men sional simulation is de- 
scribed by lYorke & Bodenheimerl d!999l) . The simulation con- 
sists of multiple nested grids in order to sample areas with higher 
spatial resolution, e.g., the inner disk. The smallest scales cor- 
respond to ~ 7 AU. The central cell is treated as a sink cell 
into w hich mass and angular momentum can flow dBoss & Black 
1 1 982b . The Poisson equation for the gravi t ationa l potential is 
solved according to iBlack & Bodenheimerl (119751) . For the hy- 
drody namics and radiation transport, the scheme from Roz vczkal 
( 19851) is tak en. Artificial viscosity is included for shocks. The 
scheme from Laug hlin & Bod enheimerl ( 1 19941) is used for mod- 
eling angular momentum transport. 

The initial conditions we ad opt correspond to the J-model 
of lYorke & Bod enheimerl (1 19991) : a 1 M Q isothermal sphere (a 
temperature of 10 K for both the gas and dust) with an initial 
radius of 6667 AU and a power law density profile with a slope 
of -2. The isothermal sphere was given a solid body rotational 
perturbation of 10~ 13 s _1 . For this set of initial conditions the 
free-fall time (tg) is 1 x 10 5 yr. The simulation runs from t — 
to t = 2.5 tff. 

2.2. Chemical model 
2.2.1. Gas phase chemistry 

For computing the time-evolution of abundances of various 
species we make use of the UMIST 2 07 da tabase^] of gas 
phase chemical reactions, Wood alTet alJ (I2007I) . This updated 
network consist of 420 species linked by 4573 reactions. We use 
the RATE06 dipole version of the network which has dipole- 
enhanced ion-neutral rate coefficients (applicable at low temper- 
atures). The initial abundances are shown in TableQ] 



2.2.2. Gas-grain interactions 

We extend the gas phase chemical network with gas-grain inter- 
actions: freeze-out and desor ption rates are calc ulated explicitly 
for each neutral species (e.g., Ber gin et alJ 1995b . Species frozen 
out create ice mantles on the dust grains. The depletion (freeze- 
out) rate, fcdep, for neutral species can be written as 



^dep 



ji a vSrigr 



(1) 



with a the mean grain radius, v = JSksT gas I 'nm(X) the mean 
thermal velocity of the gas, S the sticking coefficient, m(X) the 
mass of species X, and n gr the grain number density. The stick- 
ing coefficient S is a function of r gas and Tdust, depends on the 
species accreting, the velocity of the gas, the interaction energy 
between the grain surface and the gas species, and the excita- 
tion of the phonon spectrum of the grain. We ignore most of 



these complications and adopt a sticking probability of unity, ex- 
cept for atomic h ydrogen for w hich S is described by a function 
of r gas and rd ust (lTielensll2005l) . Past studies hav e also consid- 
ered l ower values for the sticking coefficients (e.g., Aikaw aet al.l 
1996). Lower limits for the sticking coe fficients for CO, O2, an d 
N2 at Tdust — 15 K w ere measured by Bissch op et al.l ((2006); 
lAcharyya et al.l (120071) . The reported values were 0.9, 0.85, and 
0.87 respectively, with a typical error of 0.05, so sticking co- 
efficients which are significantly lower than unity do not seem 
to be justified, except at higher temperatures of T > 10 2 K. 
Temperatures in our model are generally lower than this value 
so we choose not to lower the sticking coefficients. We take a 
"standard" mean grain radius of 0. 1 fim, and a grain abundance 
of 1.33 x 10~ 12 relative to n H = «(H) + 2n(H 2 ), with n(H) the 
number density of atomic hydrogen, and n(H2) the number den- 
sity of molecular hydrogen. 

When a species (X) freezes out it becomes bound. The val- 
ues for the binding energy, E b (X), onto the surface are uncer- 
tain. Just a few have be en measured in lab oratories, examples 
are CO, N2, O2 (see also Collings et al. 2004). The binding ener- 
gies depend on the type of grain mantle (e.g„ bare Si02, CO ice, 
H2O ice, etc.). To complicate things even further, grai n mantles 
are mi xtures of ices and evolve with time. Following iLee et all 
(2004), binding energies of a bare SiQ2 mantle are used as a 
starting point (Hasegaw a et alJ[T992t IHasegawa & Herbst 1993, 
and references therein). The binding energies onto an H2O man- 
tle (polar) and onto a CO mantle (a-polar) are greater and smaller 
than that of a Si02 mantle by factors of 1.47 and 0.82, respec- 
tively. Note that these ratios are just rough estimates and could 
vary between species. Here we ad opt CO mantle binding ener- 
gies following Aikawa et al. ( 1996), i.e., we multiply the binding 
energies by a factor of 0.82. For some species no binding ener- 
gies were found in the literature. For species which are chemi- 
cally similar to species with binding energies given, £), is inter- 
polated. Other wise E\, is estimat ed as £\, = 50A[, with A{ the 
atomic weight dWiebe et al.l l2003) and E b in units of Kelvin. 

Species can be removed from grain mantles by various 
mechanisms. It is not completely clear yet which mechanisms 
operate and under what conditions, but at least thermal desorp- 
tion (sublimation) is thought to be the most important one in pro- 
tostellar environments. The thermal desorption rate, £ su bim, can 
be written as 



ublm 



v(X) exp 



-E b (X) \ 



(2) 



with v(X) the vibrational frequency of X in its binding site, 
and Tdust the dust temperature. The vibrational frequency is the 
characteristic timescale for a species to acquire sufficient energy 
through thermal fluctuations in order to sublimate. With the har- 
monic oscillator approximation, the vibrational frequency of X 
in its binding site is given by 



v(X) 



2n s E b (X) 
"y n 2 m(X) 



(3) 



http://www.udfa.net 



with n s is the surface density of adsorption sites on a grain. The 
value of n s is estimated to be ~ 1 .5 x 10 15 cm -2 (IHasegawa et al.l 
119921) . 

A second desorption mechanism which is taken into account 
is cosmic r ay induced desorption. Th is mechanism was first 
proposed bv lWatson & Salpeterl (|1972|) . Energetic nuclei might 
eject molecules from grain surfaces by either raising the temper- 
ature of the entire grain or by spot heating near the impact site. 
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Cosmic ray heating is dominated by heavy ions because they 
deposit much more energy compared to lighter more abundant 
energetic nuclei. It is important in cold molecular cores as the 
temperature s there are too low for therm al desorption. The for- 
mulation o f lHasegawa & Herbstl dl993l) is used to calculate the 
cosmic ray desorption rate 



Table 1. The initial abundances 



3.16- l(T 19 v(X)exp 



-E b (X) 
k B ■ 70 [K] 



(4) 



2.2.3. Dissociative recombination on grain surfaces 



Under normal conditions ins ide a molecular cloud, most g rains 
will be negatively charged (Umebavashi & Nakanol l 1980). It is 
assumed that when a positive ion collides with a grain it will 
undergo dissociative recombination, with the same branching 
ratio as a dissociative recombination reaction in the gas phase. 
The products of these reactions return to the gas phase imme- 
diately. After the positive ion and the grain surface have been 
neutralized, an electron will stick immediately to the grain sur- 
face restoring its negative charge. If we assume that all grains 
have a negative charge then positive ions will collide more often 
with grains than their neutral counterparts. For single ionized 
species the collision rate is enhanced by a factor C; on . The rate 
coefficient for a dissociative recombination on a grain surface is 
then written as 



kdr = ai7ra 2 vn gv C k 



with 



Cion =51 + 



1.671 x 10" 



gas 



(5) 



(6) 



and a and T^s in cgs units dUmebavashi & Nakanol 1980; 
iRawlings et al.l [l992l) . The sticking coefficient for all ions is 
taken unity. The summation in Eq.[5]is over all possible recom- 
bination channels for a single ion; the quantity a-, denotes the 
probability for a particular channel (branching ratio). 

2.2.4. Grain surface chemistry 

Chemical reactions between species on grain surfaces can be 
of considerable importance: for example the formation of H2 
on grain mantles is the dominant mechanism to form H 2 . For 
other possible grain surface reactions the situation is less clear 
but abundances of some species are hard to explain without sur- 
face reactions as the gas phase producti on is too low to account 
for the observed abundances (see also Herbst 2005). We have 
extended the gas phase che mical network by grain s urface reac- 
tions using the approach of Hasegaw a et al.l d!992l) . Two reac- 
tants wander over the surface of a dust grain until they find each 
other and (possibly) react. This wandering/diffusion can either 
happen by thermal hopping over energy barriers or by quantum 
mechanical tunneling (only efficient for light species). The dif- 
fusion rate of species i is given by the largest of the thermal 
hopping or tunneling rate 



^"thermal — 

v(X)exp[E B (X)/k B T d 



quantum 



= v(Xy l exp[(2a/h)(2m(X)E B ) 1/2 ] 



(7) 
(8) 



For the energy barrier against diffusion, E B , w e take a value of 
30% of the binding energy, i.e, E B = 0.3£b (Hasegaw a et al.l 
|1992|) . The reaction rate between two species i and j is given by 



Species 


Abundance" 




(Relative to the number of hydrogen nuclei) 


H 


0.0 


H 2 


0.5 


He 


9.75(-2) 


O 


1.80(-4) 


C + 


7.86(-5) 


N 


2.47(-5) 


Mg + 


1.09(-8) 


Na + 


2.25(-9) 


re 


O "7/1 i fl\ 

2. /4(-y; 


S + 


9.14(-8) 


Cl + 


1.00(-9) 


p + 


2.16(-10) 


F 


2.00(-10) 


Si + 


9.74(-9) 


e~ 


7.87(-5) 


Grains 


1.33C-12) 



" abundances are given as a(b) = a x 10* 

The in itial abundances, most of them were t aken from lAikawa et al.1 
( 120011) . the grain abundance was taken from lHasegawa et alj dl992h . 
For the Fluorine a bundance we take a fifth of Chlorine abundance 
dZiurvs et all 19941) . 



with N s the number of surface sites on a grain (~ 10 ), <cy = 
exp[-2(a t /h)(2jjE a ) 1 ^ 2 ], with a t the tunneling barrier (1 A), yu 
the reduced mass, and E a the activation energy for the reac- 
tion. For most reactions E a — so that Kg = 1. The se t 
of surface reactions was taken from Hasega wa et al.1 d!992l) ; 
lHasegawa & Herbstl d 1993b . model N(4200K). This set of reac- 
tions was complimented by that of ICaselli et al.l d2002l) . but in- 
cluding only the non-deuterated species. Activation energies for 
some specific reactions taken were: 2500 K for H + CO — > HCO, 
2500 K f or H + H 2 CO -> C H 3 Q, and 176 K for OH + CO -> 
C O7 + H dChang et al.l[2007h . For CO + O -» C0 2 we took 1000 



K (d'Hendecourt et al. 



1985). 



Th e sur face rates were alter ed according to Stantc heva et al.l 
d200ll) and ICaselli etail d2002l) in order to cope with the dis- 
crete aspects of grain surface chemistry. This so called "modified 
rate approach " has been tested aga i nst a m ore detailed stochas- 
tic method bv lStantcheva & H erbst (2004) and shows reasonable 
agreement with it. In short, we compare the reaction rate coef- 
ficient (Eq.|9]) with that of the accretion and desorption rates. If 
the reaction rate is the smallest amongst the three the rate is not 
altered. In all other cases, the diffusion rate is replaced by the 
larger of the accretion and desorption rates. The net effect is to 
slow down the reactions which include at least one light species 
(i.e., atomic hydrogen). While not entirely correct, stochastic 
methods are still not available for large chemical networks and it 
is at least preferable the use the "modified rate approach" instead 
of the ordinary surface reaction rates. 



3. Method 

We integrate our initial abundances (Table [T) for 10 5 yr us- 



(9) 



ing typical conditions for a dense interstellar medium (ISM): 
n — 10 4 cm -3 , J gas = Tdust = 10 K, and an extinction corre- 
sponding to Ay = 5. For this integration time, the gas phase 
abundances and grain mantle composition show a good match 
with observed abundances and ice-compositions in dark molec- 
ular clouds. The updated abundances then serve as starting abun- 
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Fig. 1. Distribution of trace particles at t — 2.5 fff. The parti- 
cle density increases inwards to sample these regions at a higher 
spatial resolution. In this example, the particles are color-coded 
according to the fractional CO abundance. For each of these 
trace particles we solve for the chemical abundances using Eq. 
[TUl and QT| The trace particles are more or less regularly dis- 
tributed as the initial starting positions, {Rq,zq) at t = 0, were 
chosen such that at t = 2.5 % they had moved to their respective 
location in order to properly sample the abundances throughout 
the core. 
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Fig. 2. Distribution of CO throughout the core at four differ- 
ent times during the simulation for Model S. The black contour 
shows the outline of the disk. The disk is defined as the region 

where -yjv\ + v 2 R < v^, with the rotational velocity. 



face fractional abundance x s (j) for species i 

Zj K U x (^ x U) n a 2 ~ x(i) 2^ Kijx( j)n H2 



dx(i) 
dt 



(10) 



' j j 

-k dep (i)x(i) + [£ su bim(0 + k a -(i)]x s (i) , 



dances for the chemical simulation of the collapsing core. Before 
the actual collapse we also take a static core configuration into 
account which lasts 3 x 10 5 yr (i.e., one dynamical timescale and 
roughly the lifetime of a prestellar core). The abundances are 
then further integrated corresponding to this static configuration, 
but now with the temperature and density given by the first snap- 
shot (f = 0) of the hydrodynamical sim ulation. This approach 
is similar to that of lAikawa et al.l (2008). For direct cosmic -ray 
ionization reactions and cosmic-ray-induced photo-reactions we 
use the UMIST 2007 ionization rate ( of 1.4 x 10~ 17 s" 1 for 
H2. For the photodissociation reactions we also use the rates as 
given by the UMIST 2007 database, corresponding to a stan- 
dard unshielded interstellar ultraviolet radiation field of about 
10 8 photons cirT 2 s _1 between 6 and 13.6 eV. For the extinction 
we take a constant value of Ay = 15 throughout the core during 
the whole simulation which implies that the UV-radiation field 
has little effect on the chemistry. The cosmic ray ionization rate 
is left unaltered. 

The output of the hydrodynamical simulation serves as an 
input for the physical conditions needed in the chemistry. As in 
Paper I, we populate the computational domain of the first snap- 
shot (f = 0) with trace particles, see Fig. [TJ The particles have 
no interaction with each other or the medium. The dynamics of 
these trace particles is discussed in Paper I. Instead of solving for 
the CO abundance only (using a very simple scheme), we com- 
pute abundances using a full set of chemical reactions (including 
gas-grain interactions and surface reactions). 

For each trace particle we solve the following differential 
equations for the gas phase fractional abundance x(i), and sur- 



dx s (i ) V - ! V - ! v— 1 

— — = 2^2^ fy/ x *(0*sO>H 2 - x s (i) 2^ kijX s (j)n Hl (11) 

1 j j 

+k dep (i)x(i) - [£ su bim(0 + k CI (i)]x s (i) , 

with Kij the gas phase reaction rates between species / and j and 
kjj the grain surface reaction rate between species i and j. In both 
equations the first and second terms correspond to the gas phase 
formation and destruction of the species respectively. The last 
two terms correspond to the depletion and desorption processes. 
Since the particles move with the flow of the gas the absolute 
abundances evolve according to the molecular hydrogen density 
nn 2 oc n(R, z, t). The set of differential equation s is solved using 
DLSODE dRadhakrishnan & Hindmarshlll993l) . Thus, for a sin- 
gle trace particle the calculation goes as follows: a trace particle 
starts at a position (Rq, zo) within the core with the abundances 
found from our initial calculations. Rate coefficients (Kij and kij) 
are calculated using the density and temperature at that position. 
The abundances in this trace particle (fluid element) are updated 
by integrating the differential equations of our chemical network 
(Eq. [10] and [TTl i. The integration time is equal to the time be- 
tween two neighboring snapshots of the hydrodynamical simu- 
lation. The position of the trace particle is then updated, with the 
velocities given by our 2-dimensional hydrodynamical model. 
This process is repeated until the end of the simulation (~ 2100 
snapshots) or when the particle falls onto the protostar. We then 
proceed with the next trace particle. 

We limit the number of particles to about 700 in order to 
carry out the computations in a reasonable amount of time. The 
trajectories of these 700 trace particles are a subset of that fol- 
lowed in Paper I. Out of 9 x 10 5 particles, 700 are chosen such 
that at particular time steps (f = 0.0, 0.5, 1.5, and 2.5 %) these 



6 



R. J. van Weeren et al,: Modeling the chemical evolution of a collapsing prestellar core 



icr 4 


t = 


tf, 








1(T 6 








c I u 






r /~~ 




HCNT- -,. ^ " 


> 10 -to 








_CN _ 




C 1Q -12 








- HCO* 


io- u 

10~ 16 


















- N 2 H + 


t=o.5 t ff : 


1 (T 4 






10~ 6 












c 10- 8 

3 io- 10 
= 1 - 12 


/ 


/ V.'' A| 

'-- / 








io- 14 

10 -16 


t= 


1.5 t ff - 




t=2.5 t„; 



0.01 



0.10 1.00 
R [x 1 3 AU] 



0.01 



0.10 1.00 
R [x 1 3 AU] 



Fig. 3. Model S. The radial abundances through the disk (R, z = 
0) at four different times during the simulation. The line styles 
corresponding to different species are indicated in the top right 
panel at t = 0.5 tff. Labeling is the same in the other three panels. 



particles are well distributed spatially throughout the core with 
respect to the varying physical conditions (i.e., the density of the 
particles increases radially inwards to properly sample the inner 
regions). We thus follow 4 different subsets of about 700 parti- 
cles each. The particles move with the flow of the gas so they 
are irregularly distributed with respect to their locations (R,z). 
Therefore a mesh of triangles is constructed with data points 
(corresponding to the positions of the particles) at the vertices 
of the triangles. The mesh of triangles defines a piecewise-planar 
interpolating function, and we can thus grid the abundances onto 
a regular 2-dimensional grid. The spatial resolution varies from 
~ 20 AU in the core center to ~ 500 AU in the outer envelope. 

In order to asses the effects of grain surface reactions on 
the abundances of species we run the same simulation again but 
without the grain surface reactions, referred to as Model G. The 
model which includes grain surface reaction we will refer to as 
Model S. 



4. Results 

4. 1 . Abundance profiles & column densities 

The 2-dimensional abundance distributions within the core for 
CO (Model S) are shown in Fig. [2] We start with a centrally 
depleted core but as the temperature starts to rise in the center 
a sublimated zone develops. The depleted zone shrinks: as the 
cloud collapses the density in the envelope decreases so that the 
freeze-out rate drops. The sublimation radius is located at ~ 500 
AU. At the later stages CO freezes out again in the disk due to 
the increase in density in this region. Qualitatively the results are 
similar to that in Paper I, but the total CO abundance of 5 x 10~ 5 
is lower than the cosmic abundance of 2 x 10~ 4 used in Paper 
I . This is due to the fact that CO ice on grain mantles effi- 
ciently transforms into other species (e.g, CO2 and CH3OH) by 
grain surface reactions, so the amount of depletion seems higher. 
Model G has a higher CO abundance (~ IO -4 , similar to that in 
Paper I), because no surface conversion of CO takes place, but 
apart from this scaling the abundance distribution looks like that 
in Fig. [2] 

For other neutral species the pattern is similar, gas-grain in- 
teractions determine to a large extent the abundance evolution of 
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Fig. 4. Model S. The radial abundances through the disk (R, z = 
0) at four different times during the simulation. Labeling of 
species as in Fig. [3] 
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species during prestellar core collapse. We start with a high level 
of depletion (several orders of magnitude locally), and when the 
temperature rises in the center a sublimated zone develops. In 
the outer parts of the envelope the density drops, resulting in the 
increase of the gas phase abundances of species with a relatively 
low binding energy. Within the disk, species with a low binding 
energy freeze-out again. Species with a higher binding energy 
(> 1200 K) are able to enter the disk without sublimating. Gas 
accretes onto the disk around ~ 500 AU (at the centrifugal ra- 
dius) well outside the sublimation fronts for these species. Some 
ices (e.g., H2O and CH3OH) in the disk are thus (partly, as they 
are also being produced within the disk) of primordial origin. 
The abundances of charged species are determined by the bal- 
ance between destruction and formation. Parent species, and/or 
destructors often have gas-grain interactions so the abundances 
of charged species are influenced by gas-grain interactions in- 
directly. The 2-dimensional abundance distributions for various 
species are included in the online Appendices lAl (Model S) and 
EI (Model G). 

Abundances of various species as function of radius (R, z = 
0) are shown in Fig. [3] and [4] Abundances perpendicular to the 
disk plane (R = Q,z) are shown in Fig.[5]and[6] We have chosen 
to show only the results for Model S as we will discuss the differ- 
ences with Model G in the next sections (with respect to global 
envelope abundances). The radial profiles, cutting through the 
disk midplane, are difficult to compare with previous work due to 
the fact that the one-dimensional models do not include rotation 
and hence do not form a disk. We can compare the a bundance 
profile s pe rpendicular to the di sk plane with that of iLee et~aT] 
(2004J) and lAikawa et alj d2008l) . The underlying physical mod- 
els and chemical networks differ so that a one-to-one compari- 
son of the results is not possible, instead we will try to identify 
differences or similar patterns in the abundance profiles. 

We can see various abundance fluctuations for species espe- 
cially at smaller radii (< 1000 AU) or near their sublimation 
front. In particular, species often show changes in the abun- 
dances near the CO front, indicating that their formation or 
destructio n is dependent on the CO abundance, as was also 
noted by ILee et all (120041) . For example the HCO + abundance 
increases near the CO sublimation front, while the N2FF abun- 
dance slightly decreases. Further inwards the abundances of 
these two spe cies change again du e to other gas phase reactions. 
Compared to lAikawa et al] (12008). the NH3 abundance shows a 
smaller increase at its sublimation front, about one order of mag- 
nitude compared to three orders of magnitude. This is caused by 
our lower binding energy of 887 K compared to 5534 K. We 
have chosen the lower value for the binding energy as we as- 
sume that the dust grains have an a-polar CO-mantle. This value 
may be correct at low temperatures (T < 20 K), but could be 
higher in the inner regions as most of the CO mantle has subli- 
mated and we are left with a mantle consisting mostly of H2O 
ice. The CS abundan ces are two orders of magnitude lower then 
lAikawa et"aT1 d2008l) . while for HCN they are higher by about 
the same amount. The HCO + , N2H + , and H2CO abundances 
are similar. Note that when we compare the results for Model 
G there is better agreement for some species, in particular for 
HCN and CS. The difference between Model S and G will be 
discussed in the next sectio ns. 

Compared to lLee et alj d2004h . CS, HCN, H 2 CO, and NH 3 
abundance fluctuations are smaller, in par ticular at small radi i 
and near sublimation fronts. As argued by Aika wa et al] ((2008), 
these large fluctuations could arise from the fact that the number 
of species (~ 80) and chemical reactions (~ 800) are signifi- 
cantly lower compared to our network. In a small reaction net- 
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Fig. 7. Column densities for various gas phase species at four 
different times during the simulation, t — (solid lines), t = 0.5 
tfs (dotted), t — 1.5 % (dashed), and t = 2.5 % (dash-dotted). 



work, the sudden change of the abundance of one species may 
easily change the abundance of other species. In larger networks, 
there are more species included in the formation and destruc- 
tion channels of a particular species so a sudden change in the 
abundance of one species does not necessarily propagate towards 
other species. 

We have calculated vertically integrated column densities for 
various species, (iVx = J nx(z)dz), which are shown in Fig. [7] 
At t = 0, these species show flat profiles within ~ 1000 AU. 
After the collapse neutral species show an increase in the column 
density at the sublimation front. For HCO + , there is an increase 
around 300 AU (at t — 1 .5 and 2.5 fff) but further inwards the col- 
umn density decreases again and remains flat. For N2H + the sit- 
uation is similar. CN shows relatively large variations in its col- 
umn density inward s of ~ 200 A U. The column densities from 
lAikawa et all (l200ll 120031 120051) can be compared to our first 
snapshot (f = 0). The results are comparable for most species, 
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flat profiles towards the core center and slightly dropping for 
R > 5000 AU. Detailed comparison is not possible as our un- 
derlying physical model is rather different. Column densities 
in the disk increase ab o ut three to four orders o f magnitude as 



in 



Aikaw a et all dl996h : lAikawa & Her bst (2001); Aikaw a et all 



(2002) but are higher overall due to our massive disk of 0.4M o . 



4.2. Comparison with observations of envelope abundances 



J0rge nsen et al.ld200 2, 2004, 2005a) presented a sample of abun- 
dances of various species in the envelopes of Class and Class 
I objects, and in a few prestellar cores. Abundances were de- 
rived using observed line intensities and Monte Carlo radiative 
transfer modeling. In paper I, our CO abundances show the same 
trend as J0rge nsen et al.l d2002l) . i.e., a rising CO abundance for a 
decreasing envelope mass. The global envelope abundance of a 
species in the simulation excludes the disk. The envelope is de- 
fined as the region where ^jvl + v 2 R > with v$ the rotational 

velocity. The abundance is calculated by dividing the total num- 
ber of molecules of a species in the envelope by the total number 
H2 molecules in the envelope. We will now discuss the results for 
CO and the oth er species which were present in the sample of 
J0rg ensen et all We have also included the results from Model 
G, which does not include grain surface reactions. By comparing 
the results of Model G and S we can determine the influence of 
surface reactions on the abundances. 

4.2.1. CO 

The modeled CO abundances as function of envelope mass are 
shown in Fig. [8] The abundance increases from 7 x 10~ 6 to 
4 x 10~ 5 , for Model S. Overall the CO abundance is lower by 
about a factor of 10 compared to Model 3 in Paper I (includ- 
ing cosmic ray desorption and with E\, = 960 K). Some of the 
Carbon (and Oxygen) ends up in other species than CO. In par- 
ticular, CO on grain surfaces is transformed efficiently towards 
CO2, H2CO, CH3OH and other species. If the grain surface reac- 
tions are slowed down, we will get a better match to the observa- 
tions. In fact, CO2 ice and H2CO (gas) abundances are slightly 
overproduced in our simulation. The results of the simulation 
without grain surface reactions, Model G, are very similar to that 
in Paper I. This confirms that the CO abundance is mostly deter- 
mined by gas-grain interactions and grain surface reactions, and 
not by the gas phase chemistry during prestellar core collapse. 



4.2.2. HCO + 

The main formatio n and destruction m echanisms for HCO + in 
the envelope are dJ0rgensen et al.l l2004) 



Ht + CO -» HCO + + H 2 



HCO + +e + -^CO + H 



(12) 



(13) 



Before collapse, the HCO + abundance decreases radially in- 
wards. This is caused by freeze-out of its parent molecules 
(mainly CO) and by dissociative recombination of HCO + (on 
grain surfaces). Within the sublimation front HCO + shows a de- 
crease in its abundance at smaller radii because of a lower 
abundance. In Fig. [9] we have plotted X Hrn + a gainst envelope 
mass for 18 objects from J0rgens e~n et al.l d2004). Our evolution 
of the HCO + abundances agrees quite well with the observations, 
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Fig. 8. Global envelope abundances of CO versus envelope mass. 
The bullets "•" connected by solid lines are for Model S, 
which includes grain surface reactions. The bullets connected 
by dashed lines are for Model G without grain surface chem istry. 
Data for various objects from J0rgens e~n et al.l (|2002, 2004]) are 
overplotted. Class I objects are indicated by "0", Class objects 
by and prestellar cores by 
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Fig. 9. Global envelope abundances of HCO + versus enve lope 
mass. Data for various objects from J0rge nsen et al.l d2004l) are 
overplotted. Symbols and line styles are defined in Fig. [8] 

increasing from 10~ 9 to 10~ 8 , indicating its chemistry is rela- 
tively well understood. The model without grain surface reac- 
tions shows a somewhat steeper evolution but overall is very sim- 
ilar, which also shows that the HCO + abundance is not greatly 
influenced by grain surface reactions. 

4.2.3. N 2 H + 



The main formation mechanism for N2I-F ( J0r gensen et al.l 
I2004 is: 



H+ + N 2 -> N 2 H + + H 2 . 

The main destruction mechanisms are: 

N 2 H + + e " -» N 2 + H 

N 2 H + + CO -> N 2 + HCO + . 



(14) 

(15) 
(16) 



The N2I-F abundances decrease radially inwards before col- 
lapse. This is caused by freeze-out of its parent molecules 
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Fig. 10. Global envelope abundance s of N2H + versus env elope 
mass. Data for various objects from J0rge nsen et al.l (120041) are 
overplotted. Symbols and line styles are defined in Fig. [8] 
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Fig. 11. Global envelope abundanc es of CS versu s env elope 
mass. Data for various objects from J0rgensen et al. (2004]) are 
overplotted. Symbols and line styles are defined in Fig. [8] 



(mainly N2) and dissociative recombination of N2H + (on grain 
surfaces). At t — 0, the N2I-F abundance drops slower inwards 
than the CO abundance. Within the CO sublimation front, N2H + 
is destroyed by CO, at later stages this happens throughout the 
envelope. The low H^ abundance in the center also slows down 
N2FF production. However, in the envelope the net abundance 
slightly rises due to a decr easing density. Mod el G shows lower 
abundances (also noted by Aikawa et al. 2003) and steeper evo- 
lution, both in disagreement with the observations. The higher 
abundances of about a factor two in the model with surface re- 
actions are the result of the lower CO abundance in combination 
with the effects of grain surface formation of N2 and subsequent 
desorption of this parent molecule to the gas phase. We conclude 
that the model with grain surface reactions gives better results for 
N 2 H + . 

4.2.4. CS and SO 

The sulfur-bearing species CS and SO show relatively low abun- 
dances of 10" 11 - 10" 10 during the simulation. These molecules 
are heavily depleted before collapse due to relatively high bind- 
ing energies of 1780 K for both species. The abundances of 
both species increase by about an order of magnitude during the 
collapse due to the lower density (outer envelope) and higher 
temperature (center). For SO, in the center a spherical central 
sublimated zone develops (f = 0.5 ts) with an abundance of 
~ 4 x 10" 10 . Above the midplane of the disk at a radius of ~ 200 
AU the abundance rises towards ~ 3 x 10" 9 . For CS the pattern 
is similar, an abundance of 4 x 10" 11 in the center at t — 0.5 
tg, the abundance peaks around t — 1.5 tg at 1 x 10" 8 , above 
and within the disk at 100 AU. The enhancement is however too 
localized to have a major influence on the global envelope abun- 
dances. Observed global abundances are more then an order of 
magnitude higher and the abundances for Model G agree with 
the observations. This difference could be caused by grain sur- 
faces reactions converting both species towards other species. If 
we would lower the rates for these reactions surface abundances 
of CS and SO increase, and as a result gas phase abundances 
should increase as well. The increase or decrease in the abun- 
dances of other species, linked to reaction networks of CS and 
SO, could also have an effect. Abundances of both species can 
also be increased by the impact of outflows which processes the 



gas (e.g., Bach iller & Perez GutierrezTT9 97). As outflows are not 
included we could underestimate the abundances. 

It has been suggested t hat sulfur-bearin g species might pro- 
vide a chemical clock (e.g., Ruffle et al. 1993) • In the case of low 
mass p rotostars this was further investigated bv lBuckle & Fuller! 
d2003l) . The chemical clock works as follows: large amounts of 
H2S are released from grain mantles by the rising temperature 
in the inner envelope. H2S is converted via several reactions 
steps into SO and SO2 on timescales of 10 4 - 10 5 years. The 
abundances of these species drop at later times as they are con- 
verted into CS, H 2 CS, and OCS. Model S shows such a trend, 
at later times the SO abundance slightly drops and the CS abun- 
dance continues to increase. By comparing the abundances of 
species like CS, SO, SO2, H2S, and H2CS the timescale since 
the rel ease of H2S can hopefully be constrained. lBuckle & Fuller! 
(2003) found that to reproduce the observed abundances the ini- 
tial H2S abundance should be > 10" 8 . They did not compute the 
initial abundance of H2S. Instead, they estimated the abundance 
from observed gas phase and ice abundances, since no grain sur- 
face chemistry was included in their model. In our simulation 
the surface abundance of H2S is ~ 7 x 10" 8 before collapse in 
agreement with their estimates that Xh,s > 10" 8 , although the 
abun dance is on the low s ide. However, observed abundances of 
SO dJ0rgensen et al.ll2004l) do not confirm the picture described 
above, i.e., they did not find systematically higher SO abun- 
dances for Class objects compared to Class I objects. We do 
also not see such an effect in our model without grain surface re- 
actions. More observations of various sulfur-bearing species are 
needed to determine if they could be used as a chemical clock. 

4.2.5. HCN, HNC and CN 

HCN and HNC abundances observed by J0rgen sen et all (2004) 
do not show a clear correlation with envelope mass. CN abun- 
dances increase with decreasing envelope mass. This trend for 
CN is also apparent in our simulation, although abundances are 
an order of magnitude too high compared to the observations. 
CN abundances are too low by one order of magnitude for our 
model without grain surface reactions. HCN and HNC abun- 
dances are very similar and vary from 10" 8 - 10" 7 . Initially, 
abundances rise, at t > 1.5 tg the abundances remain more or 
less constant. HCN abundances are an order of magnitude too 
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Fig. 12. Global envelope abundanc es of SO versu s env elope 
mass. Data for various objects from J0rge nsen et all (120041) are 
overplotted. Symbols and line styles are defined in Fig. [8] 
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Fig. 14. Global envelope abundance s of HNC versus env elope 
mass. Data for various objects from J0rge nsen et al.l d2004l) are 
overplotted. Symbols and line styles are defined in Fig. [8] 
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Fig. 13. Global envelope abundance s of HCN versus env elope 
mass. Data for various objects from J0rge nsen et al.l (120041) are 
overplotted. Symbols and line styles are defined in Fig. [8] 
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mass. Data for various objects from J0rge nsen et al.l d2004l) are 
overplotted. Symbols and line styles are defined in Fig. [8] 



high, HNC abundances around two orders of magnitude. Both 
are produced in the gas phase but also by grain surface reactions. 

The gas phase production of HCN, HNC, and CN depends 
on the abundance of nitrogen-bearing species. Several critical 
reactions involving nitrogen-bearing species have large uncer- 
tainties in the rate coefficients, caused by the lack of theoretical 
or experi mental verification of the existence of reaction barri- 
ers ("e.g.. iFlower et al.ll2006l) . Furthermore, the elemental C:0 
abundance ratio as well as the n(N)/n(N2) ratio has an effect on 
the nitrogen-bearing chemistry as was shown by these authors. 
Another important aspect is the s t icking coefficient of atomic ni- 
trogen and oxygen. IFlower etaTJ d2006l) : lAk vilmaz e t al.l d2007l) 
found that sticking coefficients lower than unity for these species 
could explain the relatively high abundances of N2H + and NH3 
in the gas phase when other molecules, such as CO, have al- 
ready been frozen out. We conclude that the surface rates in- 
volving the production of HCN, HNC, and CN should probably 
be decreased, but that the uncertainties in the gas phase nitrogen- 
bearing chemistry may also have caused our abundances to differ 
from the observed values. 



4.2.6. HC 3 N 

Our modeled HC3N abundances provide a relatively good match 
to observed values, rising from 10~ 10 to 10 9 . Without sur- 
face reactions, abundances drop by about two orders of mag- 
nitude. Some chemical models predict an increase in the abun- 
dance when the amount of depletion of other gas phase species 
(like CO) is high ( e.g., lHiraharaetal1ll992t iRufHe et al.lll997l: 
ICaselli et al.l fl998). Observed abundances do not show such a 
trend as does our model, in fact our model predicts slightly 
higher abundances with decreasing envelope mass due to subli- 
mation of HC3N from grain mantles and lower freeze-out rates. 
J0rge nsen et al.l (120041) mention the possibility (amongst oth- 
ers) that the UV (ultraviolet) radiation from the central protostar 
could also lead to such a trend. Our results do not support this 
explanation as we have not included the UV radiation from the 
central star. In fact, HC3N produced on grain surfaces and sub- 
sequently released to the gas phase could produce the observed 
values. 

4.2.7. H 2 CO and CH 3 OH 

Observed H2CO abundances (J0rgensen et al. 2005a) show a 
clear anti-correlation with the envelope mass. H2CO freezes out 
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Fig. 16. Global envelope abundance s of HC3N versus env elope 
mass. Data for various objects from J0rge nse"n et al.l (120041) are 
overplotted. Symbols and line styles are defined in Fig. [8] 



at high densities and low temperatures. Sublimation takes places 
as the temperature rises in the inner regions and the freeze-out 
rate drops as the density in the envelope decreases. Our mod- 
eled abundances confirm this picture. The model without grain 
surface reactions produces abundances about an order of magni- 
tude lower. 

Our modeled CH3OH abundances show large variations, in 
the prestellar core observed gas phase abundances are very low 
(10~ 12 - 1CT 11 ) because of the large binding energy. Due to 
the rising temperature and subsequent sublimation, the CH3OH 
abundance quickly increases in the center, at later times the 
abundance drops again. The model without surface reactions 
fails to reproduce the observed abundances by about three to 
four orders of magnitude. Most of the CH3OH is produced by 
grain surface reactions, where CO is converted via several steps 
towards CH3OH which then sublimates into the gas phase. It 
is also believed that CH3OH abundances and to a lesser extend 
H2CO are infl uenced by other processes such as the impact of 
outflows (e.g. jBachilrer & Perez GutierrezHl997l). CH3OH was 
only o bserved towards a handful of objects by |j0rgensen et al] 
(2004). This makes comparison with models difficult. At least 
modeled abundances are roughly within the range of observed 
values. Increasing the number of objects with CH3OH abun- 
dances (or with tighter upper limits) would be helpful to further 
constrain the CH3OH production mechanism(s). 



4.3. Ice composition 

The ice composition of grain mantles can be derived from ob- 
servations of infrared absorption bands against a continuum. In 
dense mol ecular clouds ices are observed against background 
stars Ce.gJNummelin et alJl200lHKnez et al.ll2005l:lBergin etafl 
120051 IWhittet et alj|2007l) For Class 0, Class I, and Class II ob- 
jects the protostar itself provides the continuum against which 
we can observe ice absorption features. H2O ice is the most 
abundant ice-component followed by either CO or CO2. Other 
ices also exist but the abundances are generally lower. The H2O 
abundance on grains is a few times 10~ 4 and its abundance is 
relatively stable because the binding energy is high (we have 
adopted 4674 K). Therefore the abundances of the various ices 
are often expressed as percentage of H2O. In Model S the frac- 
tional H2O ice abundance varies between 1.55 x 10~ 4 and 1.80 
x 10 -4 outside a radius of ~ 100 AU, in agreement with ob- 



Fig. 17. Global envelope abundanc es of H2CO versus env elope 
mass. Data for various objects from J0rgens en et al] (l2005al) are 
overplotted. Symbols and line styles are defined in Fig. [8] 
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Fig. 18. Global envelope abundances o f CH3OH versus enve - 
lope mass. Data for various objects from J0rgens en et all J2004); 
ISchoier et al.l d2002h are overplotted. Symbols and line styles are 
defined in Fig. [8] Other data points indicated are the abundances 
fitted with "drop"-abundance profiles, with '+' the abundance in 
the inner region and 'x' that in the outer region of the envelope. 



served values dGibbet al.l 12000b ISchutte et al] 1 19991 and refer- 
ences therein) 

The grain mantles in Class objects have similar ice 
compositions as those in cold molecular clouds, although a 
small amount of thermal processin g seems to hav e take n 
place for a number of sources, e.g., iGerakines etail ([1999); 
Gibb et al.l d2004l) . In Class I objects evidence is found for fur- 
ther processing of ices and complex molecu l e formation e.g 
lEhrenfreund et al] (fl998l); iKeane et al] d200lb; [Alexander et al 
d2003l): I Watson et all 12004): Gi bb et all d2004l) : iBoogert et al 



( 20041); IZasowski et al.l |2007). Class II objects have heated the 
remaining envelope and outer layers of the disk above the subli- 
mation temperature of ices. Although, in some edge-on objects 
ice are ob served, probably arising f rom the cold midplane of the 
disk (e.g.. iPontoppidan et al .Il2005l) . 

The vertically integrated column densities for various ices 
(Model S) are shown in Fig. [19] For most ices the overall evo- 
lutionary pattern is similar, we start with rising column densities 
towards the center at t — in the prestellar core. Then a central 
sublimated zone develops where the column density remains rel- 
atively flat. The size of this region is determined by the binding 
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R [x 10 3 AU] R [x 10 3 AU] 

Fig. 19. Column densities of various ices at four different times 
during the simulation, t — (solid lines), t = 0.5 % (dotted), 
t = 1.5 % (dashed), and t - 2.5 ts (dash-dotted). 

energy of the species. In the cold midplane of the disk, species 
with a low binding energy freeze-out again, while species with 
a high binding energy remain frozen out. Towards the disk, col- 
umn densities increase by a few orders of magnitude for all ices. 
We will omit the results for Model G as it fails to reproduce the 
observed ice-compositions for most of our surface species, see 
also Table |2] 

We compare the modeled ice fractions at four different times 
during the simulation (t = 0, 0.5, 1.5, and 2.5 ts) to observed 
values towards the field star Elias 16 (in the vicinity of the dense 
core TMC-1), the low- mass protostar Elias 2 9, and a sample of 
Class I/II objects from Zasowski et al. (2007). The abundances 
are given relative to H2O and are averaged over the envelope 
(which does not include the disk). We have also included the 
ice composition within the disk, at a distance of 500 AU, just 
outside the CO sublimation front. Model G severely underpro- 
duces various ices like CH3OH, CO2, H2O and H2CO. This is 



explained by the fact that surface reactions are the dominant pro- 
duction mechanisms for these species. Sublimation of these ices 
then also leads to a higher gas phase abundance. Most of the 
CH4 ice is located within the disk in Model S where the density 
is high. For model G, the CH4 ice fraction is higher in the enve- 
lope. Surface reactions are thus not required to form a significant 
amount of CH4. Below we further discuss the results from Model 
S for several species. 

Modeled CO ice abundances are in rough agreement with ob- 
served values. Within the prestellar core the CO ice fraction is on 
the low side with a value of 14% while the CO2 fraction of 47% 
is quite high. Although , values of 37% and 32% were reported 
by lBoogerteTaT] (120041) towards B5 IRS 1 and HH 46 IRS. The 
CH3OH ice fraction of about 5% is in agreement with the obser- 
vations given the large scatter, the same is true for H2CO. Our 
modeled SO2 ice fraction (~ 0.02%) is on the low side, but its 
identification is uncertain. The NH3 ice fraction (~ 0.5%) is low, 
observations indicate higher fractions of about 5%, but the NH3 
identificatio ns are also somewh at uncertain. A higher binding 
energy as in lAikawa et al. (2008) would increase its fraction. As 
has been mentioned in Section l4.2.5l the large uncertainties in the 
gas phase nitrogen-bearing chemistry may also have affected the 
NH3 ice abundance. During the collapse the fraction of CO ice 
decreases from 15 to 5%, the fraction of CO2 ice also decreases, 
from 47 to 30%. The CH3OH fraction stays more or less con- 
stant while the H2CO ice fraction increases from 1 to 4%. The 
decrease in the CO ice fraction is also being observed (as dis- 
cussed above) and is the results of its low binding energy and 
the conversion into other more complex species. Slowing down 
the surface reaction CO + O — > CO2 by increasing the activa- 
tion energy of 1000 K somewhat would produce less CO2 and 
increase the CO ice fraction. An alternative would be to increase 
the barrier against diffusion (Eg) which should also slow down 
CO2 production, as well as the production of other ices. 

5. Discussion 

When we compare our gas phase abundances and ice com- 
positions to observations, it seems that the surface reaction 
rates are somewhat too high. First we slightly overproduce the 
amount of CO2 compared to CO. Our ratio of CO2 to CO is 
about three while observations often show ratios of about one 
to two. Also the H2CO gas phase abundances are a little high. 
Furthermore species like HCN, HNC are overproduced by quite 
some amount, although uncertainties in the nitrogen-bearing 
chemistry might also have influenced these results. To slow 
down the surface rates several options exist. First, the binding 
energy (E\,) of atomic hydrogen could be increased which would 
slow down surface hydrogenation reactions. We have taken a fast 
diffusion rate for atomic hydrogen (E\, = 492) so it could be that 
we have t o increase its binding e nergy to slow down the diffusion 
rate (see iRuffle & Herbstil2000h . Another option is to increase 
the barrier against diffusion Er for all surface reactions (e.g., 
£b = 0.5£b). For certain reactions it is also possible to increase 
the activation energy (E a ). Although, this will not help slowing 
down the HCN and HNC production as there are are no activa- 
tion energies involved in the production of these species but for 
the oxidation reaction of CO leading to CO2 this would be an 
option. As a first step, lowering the diffusion rate for atomic hy- 
drogen and increasing the energy barrier against diffusion will 
be the best option. 

The underproduction of CS and SO (Model S) is not com- 
pletely understood. Although CS and SO are transformed by 
surface reactions this should not have a very large effect on the 
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Table 2. Grain mantle compositions 



Species" 


Elias 16" 


Elias 29 r 


Class I/II Objects'* 


Model 


t = 0.0 


t = 0.5 t s 


t= 1.5% 


t = 2.5 tff 


midplane disk'' 


H 2 ; 


100 


100 


100 


S 


1.55(-4) 


1.65(-4) 


1.7K-4) 


1.7K-4) 


1.8(-4) 










G 


1.24(-5) 


8.29(-6) 


4.5(-6) 


4.2(-6) 


1.4(-5f 


CO 


25 


5.6 




S 


13.9 


17.0 


11.6 


5.0 


17.6 










G 


71.6 


61.1 


31.4 


13.5 


77.5 


CO, 


24 


22 


~ 12 


S 


46.9 


40.6 


31.0 


29.8 


37.7 










G 


3.6(-l) 


3.K-D 


2.2(-l) 


2.2(-l) 


4.2(-l) 


CH30H 


< 2.3 


< 4.4 


> 2-9 


S 


8.7 


7.0 


6.6 


6.7 


2.2 










G 


1.4(-5) 


2.1 (-5) 


2.9(-5) 


2.9(-5) 


7.0(-8) 


NH 3 


< 8 


< 9.2 


< 14 


S 


9.8(-l) 


3.3C-1) 


2.6(-l) 


3.4(-2) 


4.0(-l) 










G 


1.7(-2) 


l-K-2) 


2.7(-3) 


1.0(-3) 


6.0(-3) 


CH 4 


<3 


< 1.6 


~ 4 


S 


l.l(-l) 


1.7(-1) 


1.4(-1) 


1.0(-1) 


1.8 










G 


2.2 


2.0 


1.6 


1.4 


1.0 


H 2 CO 




< 1.8 


~ 2 


S 


1.1 


1.6 


3.9 


4.1 


2.9 










G 


4.5(-4) 


3.7(-4) 


4.3(-4) 


4.7(-4) 


4.8(-4) 


OCS 


< 0.2 


<0.05 




S 


4.3(-9) 


4.9(-10) 


4.6(-10) 


8.1(-10) 













G 


3.7(-3) 


4.0(-3) 


3.0(-3) 


2.7(-3) 


9.1 (-2) 


HCOOH 




<0.9 


~ 0.6 


S 


1.5 


2.6(-2) 


1.4(-2) 


1.4(-2) 


3.4(-2) 










G 


7.7(-3) 


9.3(-3) 


7.7(-3) 


7.1 (-3) 


6.9(-l) 


S0 2 






~ 0.5 


S 


2.3(-2) 


1.3(-2) 


8.7(-3) 


8.4(-3) 


1.2(-2) 










G 


7.8(-4) 


1.0(-3) 


7.0(-4) 


5.9(-4) 


8.3(-3) 



" abu ndances are given as a(b) = a x 10* 

b COlChiar et all J1995T); CQ 2 , C H 3 OH, N H 3 , and CHalKnez et alJ (l2005T) ; OCS lPalumbo et all < fl997h 

c lEhrenfreund & Schuttd feOOCl) and lBoogert & Ehrenfreund fc004l) 

rf lZasowski et all d2007l) 

e at t - 2.5 t s and a radius of ~ 500 AU 

J H 2 fractional abundance are given for the simulation: «(H 2 0)/n(H 2 ) 

g for model G we have taken an H 2 fractional abundance of 1.8 x 10 -4 for calculating the ice fraction of other species, instead of the reported 
values on this line 

Observed ice compositions as percentage of the H 2 abundance compared with model results. Every first line for a species shows the results from 
Model S (includes surface reactions) and the second line shows the results for Model G (excludes surface reactions) 



abundances as most of those surface reactions do not involve H 
atoms and are therefore quite slow. Reduced gas phase abun- 
dances of CO (and other species) could result in a different gas 
phase chemistry leading to the decrease of parent species for CS 
and SO production. More work is needed to further investigate 
this effect. 

In our simulation we have taken a constant value for the ex- 
tinction (Ay = 15) throughout the collapsing prestellar core. 
This is a good approximation during the initial phases of the 
collapse (fff < 1). Only at the outer edges of the core the results 
could change significantly due to the presence of the interstellar 
radiation field combined with lower values for the extinction. We 
are mainly interested in the central regions of the core and have 
only plotted our results for r < 5000 AU. Prestellar cores can 
also be shielded by several magnitudes of external extinction, 
for example if they are embedded within a larger dark molecular 
cloud. Our simulation can therefore be described as the chemi- 
cal modeling of an externally shielded collapsing core. At fff > 1 
the extinction in the region above the midplane of the disk drops 
below Ay < 5 so that photo-dissociation and ion-molecule reac- 
tions become important. The importance of these reactions is in- 
creased by the stronger radiation field from the central protostar. 
However, the overall effect on the global abundances of species 
is reduced as most of the mass is located within the shielded 
disk where the extinction is large (Ay > 15), even at the closest 
distance to the protostar of ~ 7 AU within the simulation. 

To include the effect of the radiation field the extinction 
should be calculated at every position within the core. From 
the density given by the hydrodynamical simulation we could 
in principle derive the amount of extinction using the relation 



^ v = 1 s9x\a 2i ~Etr? ' w ^ ^ H t ' le commn density of hydrogen 
nuclei between the protostar and the point of interest. However, 
in practice this is more complicated as a significant amount of 
extinction arises from regions close to the protostar within 7 AU, 
beyond the resolution of our simulation where we do not have 
information on the density distrib ution. The dust dest ruction ra- 
dius is located roughly at 2000 K dDuschl et alJI 19961) ~ 0.1 AU 
from the protostar which sets an upper limit on the resolution 
needed to derive the amount of extinction. By increasing the 
resolution of the hydrodynamical simulation this can be satis- 
fied. As a first approximation for the intrinsic strength of radia- 
tion field, we could take the luminosity from the protostar given 
by the hydrodynamical simulation. The protostar luminosity is 
based on the mass accretion rate plus a central core luminosity. 
Including these effects can be a starting point for future work. 
Our hydrodynamical simulation has a simplified radiative trans- 
fer, this may have resulted in deviations in the ga s and dust tem- 
peratu re compared to more realistic simulations. Jonkhei d et al.l 
(|2004|) presented a model which calculated the gas and dust tem- 
perature as well as the chemistry in the surface layers of a flaring 
protoplanetary disk. The heating rate of the gas was controlled 
by UV radiation and the gas cooling by line emission of various 
species and by collisions with dust grain. They found that the gas 
temperature exceeds the dust temperature, in the part of the disk 
that is optically thin to UV radiation, by up to several hundreds 
of Kelvin. In the optically thick part of the disk the gas tem- 
perature was found to be well coupled to the dust temperature. 
In our hydrodynamical simulation the dust isotherms ar e more 
spherical compared to the results in Jonkhei d et al.ld2004t) . Small 
changes in the dust temperature can have an important effect on 
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the sublimation rate as it is exponentially dependent on the dust 
temperature. Therefore accurate dust temperatures are important 
for calculating the position of the sublimation front if the bind- 
ing energies of the species are known. Observations of the posi- 
tion of sublimation fronts of various species could be a sensitive 
tracer of the dust temperature. The very high gas temperatures 
found above the disk also have an effe ct on the gas phase chem- 
istry. Although iJonkheid et al.l d2004l) reported that effect was 
small. However, their chemical network was rather limited and 
no gas-grain interactions were included. A temperature increase 
from 60 to 550 K in the upper disk layers (at R = 150 AU and 
z ~ 90 AU) will increase the freeze-out rate by a factor of ~ 3, 
although this may be compensated by a lower sticking coeffi- 
cient S . X-rays from the c entral protostar also h ave an effect on 
the chemistry as shown by Staub eret al.l d2005l) . Abundances of 
some species were found to be enhanced by 2 - 3 orders of mag- 
nitude. Due to smaller cross sections, X-rays penetrate further 
into the disk and envelope and therefore affect the chemistry on 
large r scales (> 1000 A U) than the UV radiation field (~ 400 
AU. IStauber et alj d2004t) ). Eventually the chemistry needs to be 
coupled to the radiative transfer (e.g., Ivan Zadelhoff et alj 120031) 
for the modeling to be self consistent and correctly address this 
issues. 

In our model we can follow the chemical composition of gas 
accreting onto the disk. Gas does not directly accrete onto the 
disk but flows around the disk upwards (see Fig. l20b before en- 
tering the disk around the centrifugal radius (see also Paper I). 
This has important implications for the chemistry as the gas is 
more exposed to UV and X-ray radiation from the central proto- 
star before entering the disk. The dust temperature at the cen- 
trifugal radius of the disk sets a sublimation temperature for 
ices accreting onto the disk. The result is that ices with a low 
binding energy (E\, < 1200 K) sublimate before entering the 
disk while ices with a high binding energy (£\, Z 1200 K) are 
able to enter the disk. The strong UV r adiation field may also 
cause photodesorption of ices. Recently, Oberg et al. (2007|) de- 
termined high yields for CO photodesorption in laboratory ex- 
periments. In our simulation at t — 2.5 tg, most of the CO ice 
(~ 75%) is not able to accrete onto the disk, as the sublimation 
front lies just outside of the centrifugal radius (at ~ 500 AU), 
while H2CO ice can enter the disk without sublimating, see Fig. 
|20l By following the chemical evolution of gas accreting onto the 
disk we can derive the initial abundances to be adopted in more 
detai led modeling of the chemistry within the disk. For exam- 
ple [Semenovetal] d2004l) took as initial abundances those com- 
puted in a dark cloud (out of which the disk eventually forms). 
However, the initial abundances are different from those in dark 
clouds as some chemical processing has already taken place be- 
fore the gas enters the disk. Using the correct initial abundances 
might improve the results of those models. 

We have chosen a standard dust grain radius of 0.1 fim 
for o ur simulation to co mpare our results with previous work 
(e.g jAikawa et al . 2008). This is an oversimplification as grains 
are expected to grow to larger sizes under the conditions in 
the simulation, especially within the disk. The main effect of 
grain growth will be a decrease in the depletion rate as the to- 
tal grain surface area decreases. The lower depletion rates also 
ha ve an effect on the n itrogen-bearing chemistry as was shown 
by iFlower et alj d2006l) . An approach mi g ht be to use the more 
sophisticated treatment by IFlower et alj (|2005), following the 
coagulation of the grains simultaneously with the chemical evo- 
lution of the medium. This is, however, beyond the scope of this 
paper. 
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Fig. 20. Left: Abundances as function of time for a single trace 
particle which accretes onto the disk. At t — 1.5% the conditions 
for the trace particle are roughly that of the dark ISM, T — 10 K 
and n = 10 4 crrT 3 . Right: Trace particle trajectory, the position 
of the particle at different times is indicated in the figure. The 
acceleration of the particle as well as the upward motion around 
the disk are visible before the particle accretes onto the disk. 



6. Conclusions 

We have modeled the chemical abundances during prestellar 
core collapse and have compared the global evolution of these 
abundances with observations. We have used trace particles, 
moving with the flow of the gas as given by a 2-dimensional 
hydrodynamics simulation, to derive the chemical evolution of 
the collapsing core. 

Gas-grain interactions determine to a large extent the abun- 
dance evolution of species during prestellar core collapse. 
Before collapse, most species are depleted due to freeze-out onto 
dust grains. When the temperature rises in the center a subli- 
mated zone develops. At the same time in the outer parts of the 
envelope the density drops, which lowers the freeze-out rate re- 
sulting in an increase in the gas phase abundance of species with 
a relatively low binding energy. Within the disk, species with a 
low binding energy freeze-out again, species with a higher bind- 
ing energy are able to enter the disk without sublimating. The 
abundances of charged species are mostly determined by the 
balance between destruction and formation. Parent species, or 
destructors often have gas-grain interactions so the abundance 
of charged species is still influenced by gas-grain interactions, 
although in an indirect way. 

We conclude that grain surface reactions are important for 
the formation of some species, in particular for CH3OH, CO2, 
H2O and H2CO. Without surface reactions most of the ice is in 
the form of CO, and the relative fraction of other ices is much 
lower than observed. An example is CH3OH where the ice frac- 
tion for Model G is lower by more than 5 orders of magnitude. 
Gas phase abundances of species likes CH3OH and H2CO are 
also too low without surface reactions when we compare global 
envelope abundances with observations. So we confirm that sur- 
face reactions are important for the production of some species. 
CH4 is produced quite efficiently in the gas phase. Only at high 
densities within the disk CH4 is produced more efficiently by 
grain surface reactions. Abundances are higher in the envelope 
if we do not include surface reactions. It could be that some CH4 
is transformed in Model S towards CH3 and C2H2 by surface 
reactions. CS and SO abundances show a better match to obser- 
vations if no surface reactions are included. 

Our method provides important information on the initial 
abundances of species to be adopted in chemical modeling of 
protoplanetary disks. The initial abundances are often derived 
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from chemical models of dark molecular clouds. However, as 
we have shown, significant chemical processing of the gas and 
the grain mantles occurs before accretion onto the disk. For ex- 
ample ices with a low binding energy (< 1200 K) are not able to 
enter the disk, while those with a higher binding energy are able 
to directly accrete onto the disk. Also gas accreting onto the disk 
(coming from z < 1000 AU) may be temporally exposed to the 
strong radiation field from the central protostar, especially at the 
later stages of the collapse when t > 2 tg. This is caused by the 
upward motion of gas around the disk, representing a bow shock 
feature. The gas then accretes above the midplane, roughly at the 
position of the centrifugal radius onto the disk. 

The method we have presented for modeling the chemi- 
cal abundances is generally applicable to other hydrodynamical 
simulations. Although, modeling the chemical evolution is com- 
putationally expensive it can be parallelized easily as the trace 
particles have no interaction with each other. This can be used 
to extend the modeling to higher spatial resolution (from the 20 
AU we have used) and towards three dimensions. This will allow 
us to probe the chemical evolution in the planet forming zone in 
greater detail. Including the effects of the UV radiation and X- 
rays is important in the later stages (fff < 1) of the collapse, 
especially for the upper disk layers and the envelope. 

Finally, modeling line profiles with RATRAN and compar- 
ing them to observations may provide more information on the 
physical conditions during prestellar core collapse. 

From an observational point of view, increasing the sample 
of objects with modeled envelope masses and abundances, will 
help to test chemical models. In particular, more CH3OH obser- 
vations as well as those from sulfur bearing species are needed. 
This will also help to find species which could be used as chem- 
ical clocks. Second, laboratory experiments of grain surface re- 
actions may further help to constrain the surface reaction and 
diffusion rates. 
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Appendix A: 2-Dimensional abundances including 
grain surface chemistry 
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Fig. A.3. Same as for Fig.|AT]but for CS 
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Fig. A.l. Distribution of HCO + throughout the core at four dif- 
ferent times during the simulation for Model S. The black con- 
tour shows the outline of the disk. The disk is defined as the 

region where -^v^ + v 2 R < v$, with v$ the rotational velocity. 
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Fig. A.4. Same as for Fig. lA.ll but for SO 
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Fig. A.2. Same as for Fig. lA.ll but for N 2 H 
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Fig. A.5. Same as for Fig.|AT]but for HCN 
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Fig. A.6. Same as for Fig.lAJlbut for HNC 
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Fig. A.9. Same as for Fig.gj]but for H 2 CO 
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Fig. A.7. Same as for Fig.|ATT]but for CN 
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Fig. A.10. Same as for Fig.|AJ]but for CH 3 OH 
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Fig. A.8. Same as for Fig.|ATT]but for HC 3 N 
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Appendix B: 2-Dimensional abundances without 
grain surface chemistry 
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Fig. B.3. Same as for Fig. IBTTlbut for N 2 H + 
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Fig. B.l. Distribution of CO throughout the core at four different 
times during the simulation for Model G. 
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Fig. B.4. Same as for Fig. EUbut for CS 
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Fig. B.9. Same as for Fig. IBTTIbut for HC 3 N 
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Fig. B.7. Same as for Fig. [BJ]but for HNC 
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Fig. B.10. Same as for Fig. |BJ]but for H 2 CO 
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Fig.B.ll. Same as for Fig.|RT]but for CH 3 OH 



